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Abstract 

Similarity renormalization group (SRG) flow equations can be used to unitarily soften nuclear 
Hamiltonians by decoupling high-energy intermediate state contributions to low-energy observables 
while maintaining the natural hierarchy of many-body forces. Analogous flow equations can be used 
to consistently evolve operators so that observables are unchanged if no approximations are made. 
I/-} The question in practice is whether the advantages of a softer Hamiltonian and less correlated 

wave functions might be offset by complications in approximating and applying other operators. 
Here we examine the properties of SRG-evolved operators, focusing in this article on applications 
to the deuteron but leading toward methods for few-body systems. We find the advantageous 
features generally carry over to other operators with additional simplifications in some cases from 
factorization of the unitary transformation operator. 
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I. INTRODUCTION 



Renormalization group methods can be used to soften interactions in nuclear systems, 
which extends the range of many computational methods and qualitatively improves their 
convergence patterns PQ. The similarity renormalization group (SRG) (2rE] does this by sys- 
tematically evolving Hamiltonians via a continuous series of unitary transformations chosen 
to decouple the high- and low-energy matrix elements of a given interaction [3, E] . In par- 
ticular, a flow equation with parameter s and generator r] s = [G s , H s ], 

= [,.,*.], a) 

unitarily evolves an initial Hamiltonian H s= o = H — T m \ + V. Choosing the flow operator 
G s specifies the SRG evolution. This equation implements the unitary transformation 

H s = U s H s=0 Uj = T rcl + V s , (2) 

which defines V s by choosing the relative kinetic energy to be invariant, and where the 
generator r) s is related to U s by 

* = ftf = -. it- (3) 

As in most previous nuclear applications, here we take G s = T rel to suppress off-diagonal 
elements of the Hamiltonian in momentum space (see SecjllJ alternative choices are discussed 
in Refs. [H] and [ID])- This decoupling leads to greatly improved convergence of the binding 
energy in few and many-body calculations [HI E]. A major advantage of the SRG relative 
to other energy-independent renormalization group (RG) methods is that the Hamiltonian 
flow equation is formulated solely in terms of the evolving Hamiltonian and does not involve 
the T-matrix, which avoids issues with solving equations in multiple channels and allows any 
convenient basis to be used; these features mean that evolving few-body forces is practical [8j 

H2]. 

To use the wave functions produced by SRG-evolved interactions to calculate other matrix 
elements of interest, we cannot in general neglect the associated change in operators. The 
evolution of any operator O = O s= o is given by the same unitary transformation used to 
evolve the Hamiltonian [HUH], 

d s = u s d s=0 ul, (4) 

which implies by differentiation with respect to s the general operator SRG equation 

^ = hA). (5) 

Although this equation can be used to find O s , it is computationally efficient to construct 
the unitary transformation directly from the eigenvectors of the evolved and unevolved 
Hamiltonian using 

u s = J2\Ms))(Mo)\ , (6) 

a 

where the sum on a is over all eigenvectors, and then to apply Eq. Q directly. In practice we 
work in a discretized basis so this sum is finite and Eq. (j4j) is a simple matrix product. In cases 
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where both methods have been used to calculate the SRG evolution, the transformations 
produced by Eqs. ^ and ^ agree up to numerical errors. However, it has been found 
that the one-step transformation produced from the eigenvectors is numerically more robust 
than the differential equation. The direct construction of U s via Eq. ^ is used to calculate 
operator evolution throughout this work. 

If implemented without approximation, unitary transformations preserve operator matrix 
elements by construction, 

(V>a(s)|6-foMs)) = (V>a(0)|O s= o|Vv(0)> , (7) 

and thus preserve the physics in the initial Hamiltonian and other operators. But do the 
advantages of the SRG evolution of Hamiltonians carry over to other operators and are there 
problems with the practical implementation of operator evolution? Equation Q implies that 
changes in the wave function are "compensated" by changes in the operator. This might 
reasonably be expected to shift around the physics while conserving the computational dif- 
ficulty, so that at best we must deal with either a simple operator and complicated wave 
functions or a complicated operator and simplified wave functions. However, the SRG evo- 
lution of the Hamiltonian generates a much simpler interaction (smoother and decoupled), 
which leads to a simpler wave function (reduced short-range correlations) [j] What about 
other operators? Could there be strong and/or fine-tuned cancellations between the evolved 
wave functions and evolved operators? Can the hierarchy of many-body contributions be 
violated? We address these questions in this article and a sequel [TO] . 

The evolution of three-body and higher-body interactions is critically important for the 
SRG and a parallel discussion is needed for other operators. To see how one-, two-, three-, 
and higher-body operators can be identified, it is useful to decompose the running SRG 
operator O s in second-quantized form. Schematically (suppressing indices and sums), 

3 S = a ] a + (Of )) aVaa + (df^) a)a ] a ] aaa + . . . , (8) 

where a) and a are creation and annihilation operators with respect to the vacuum in a 
single-particle basis. This defines (Oi 1 ^), (OT 1 ), (Of^ 1 ), ... as the one-body, two-body, 
three-body, . . . operator matrix elements in that basis at each s. The SRG evolution in 
Eqs. (jl]) and ^ is dictated by commutators involving O s and H s (which also has such an 
expansion). When they are evaluated, we see that even if initially there are only one-body 
operators, higher-body terms will appear in both O s and H s with each step in s. Thus, 
when applied in an A-body subspace, the SRG will "induce" A-body operators. However, 
we find that each (Os) is determined fully in the A = n subspace, with no dependence 
on higher-body operators. This allows us to extract the induced many-body components of 
the operator as needed [TO]. Note that for Hamiltonians with no external potentials (i.e., 
no one-body interactions), (Os) is independent of s [10J. It is also important to remember 
that input operators, O s= o, for low-energy effective theories are generally never simply one- 
or two-body operators, although these components may dominate. The question is whether 
an initial many-body hierarchy (expected from chiral effective field theory formulations of 
the nuclear interaction) is maintained by the SRG evolution. 

To avoid confusion, we note that it is also possible to normal order with respect to a 
finite-density reference state instead of the vacuum, which leads to the "in-medium SRG" 

1 Note that the interpretation of "simpler" can vary. For some Monte Carlo methods, the SRG Hamiltonians 
become more complicated in the sense of increasing non-locality in coordinate representation. 
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(see Ref. [H] and references therein). This changes the definition of the matrix elements 
and creation and annihilation operators in Eq. ^ and shifts higher-body pieces to the zero- 
body, one-body, and two-body levels. These operators are well defined but have different 
properties from those considered here; for example, even one-body operators flow. The in- 
medium SRG shows great promise as a microscopic method of deriving effective shell model 
interactions for nuclei and the study of the corresponding operators is an important topic 
for future investigation. 

In this article, we restrict our attention to the deuteron, which means only one- or two- 
body operators are relevant. Furthermore, any running with s is due to an induced two-body 
part. In practice, working in the two-body system in the center-of-mass frame makes the 
occurrence of two-body operators uneventful; there are far more consequences for A > 2. 
We defer to the sequel flO] the discussion of evolving and extracting operator components 
and embedding them in larger A systems using a harmonic oscillator basis. 

Nuclear SRG studies to date have focused primarily on the calculation of binding energies 
and phase shifts to analyze the characteristics of the SRG-evolved interactions. A limited 
analysis of the deuteron momentum distribution defined by the initial operator a^a^ has also 
been made [13J and decoupling of long-distance operators for the rms radius, quadrupole 
moment, and r _1 have been examined [7], but only for the bare operators. We expand 
upon these studies here, emphasizing the nature of SRG-evolved operators and also include 
the first calculations of deuteron charge, quadrupole, and magnetic form factors. However, 
we use these electromagnetic operators as test cases for addressing questions about opera- 
tor evolution and not yet for systematic comparison to experiment, which requires a more 
complete treatment of the s = operators. 

Other issues arise about processes with large momentum transfers, such as (e, e'p). Theo- 
retical analyses relate such experiments to nuclear momentum distributions if the impulse ap- 
proximation is assumed to be valid for a high-cutoff interaction [JS] . Calculations find nearly 
universal scaling of the high-momentum tails, which is interpreted in terms of short-range 
correlations in the nuclear wave functions. It might be thought naively that this physics is 
beyond the reach of low-momentum approaches, for which wave functions have drastically 
reduced short-range correlations. However, Eq. ^ is unequivocal: The experimental cross 
section is unchanged with SRG evolution to low momentum if no approximations are made, 
even if the evolved wave function has almost no short range correlations. But does the 
calculation become intractable because the evolution of the momentum occupation operator 
makes it too complicated (e.g., strong non-localities, too-large many-body components)? We 
begin to address this question by showing that under the relevant kinematic conditions there 
is factorization of the unitary transformation U s , which leads to significant simplifications 
and an alternative interpretation of the universal high-momentum dependence and scaling. 

The plan of the article is as follows. In Sec. |TTJ we explore whether decoupling of the 
Hamiltonian for two-body systems is mirrored in the operator flow, focusing on evolution 



of the momentum distribution as a characteristic example. In Sec. Ill we consider the 
evolution of other operators, including electromagnetic form factors in the deuteron, which 
are all found to flow to smooth, low-momentum forms. Variational calculations of the 



evolved deuteron binding energy and operator matrix elements are explored in Sec. |IV| The 
factorization of the unitary transformation operator under certain conditions is demonstrated 
in Sec. [V] along with a model calculation of the momentum distribution for A > 2. Our 
conclusions are summarized in Sec. [V]] 
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II. DECOUPLING AND THE DEUTERON MOMENTUM DISTRIBUTION 



A. Potentials and Decoupling 



In this section, we specialize the SRG evolution to a two-particle partial-wave momentum 
basis with flow operator G s = T re \. The flow equation, 



-£ = [[T t *,H.],H.] = -£ 



(9) 



(recall that T re i is chosen to remain constant) for nucleon-nucleon (NN) potentials is pro- 
jected onto each channel using 1 = | J °° q 2 dq \ q) (q\ with h — M — 1, yielding 



dV s {k,k') 
ds 



-{k 2 - k' 2 ) 2 V s (k,k') 

+ - / q 2 dq (k 2 + k' 2 - 2q 2 )V s (k, q)V s (q, k') . 
^ Jo 



(10) 



This equation is implemented in a discretized Gaussian quadrature basis as a set of coupled 
differential equations for the matrix elements (angular momentum indices in the coupled 
channels have been suppressed) that are solved numerically. 
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FIG. 1. (Color online) SRG evolution of the momentum-space 3 Si potential starting with (a) 
Argonne v 18 (AV18) from A = 15 fm" 1 to A = 1.5 fm" 1 [16J and (b) N 3 LO (500 MeV) from A = 
6ml- 1 to A = 1.5 fm" 1 [T7]. 

The flow is visualized in Fig. [T] for two representative initial potentials, Argonne v is 
(AV18) [16] and a chiral N 3 LO (500 MeV) effective field theory (EFT) potential from 
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Ref. [IT], which are currently the most commonly used NN interactions for microscopic 
nuclear structure calculations. The figures show snapshots of the potential matrix V s (k, k') 
in the 3 Si partial wave at several values of A, which is a useful alternative parameter to 
characterize the flow. It is related to s by A = s -1 / 4 . The scale at the right indicates 
the strength of the potential (in fm); Gaussian mesh weights are not shown. One can see 
that the large matrix elements between low and high momentum at the initial A's shown 
(particularly for AV18) are suppressed by A = 1.5 fm - ; that is, we see decoupling for both 
potentials. 

The decoupling seen in Fig.[T]is readily understood from the SRG flow equations. Because 



of the dominance of the kinetic energy, Eq. (10) for sufficiently off-diagonal k and k' is given 
to good approximation by 



dV s (k,k') 
ds 

which, when solved, predicts 



-(k 2 -k' 2 ) 2 V s (k,k'), (11) 



V s (k, k') « e- s{k2 - k '^ V s=0 (k, k') = e" 1 *^^ V x=ao (k, k!) . (12) 

Thus, using this generator we can see that the far off-diagonal elements of the potential 
matrix are suppressed exponentially with an approximate width given by the flow parameter 
A [6]. (If the potential is plotted as a function of k 2 , the width of the partially diagonalized 
potential is clearly seen to be well approximated by A 2 |1J.) While it is not evident from the 
figure, the potential is also very smooth. 

It is not immediately clear, however, that this decoupling will be advantageous for the 
calculation of observables other than the binding energy. If we project the operator flow 
equation 

d ^ = [[T rch H s ],d s ], (13) 
onto the partial wave momentum basis, we find 



dO s (k,k') 



ds 7T 



poo 

/ q 2 dq [(k 2 - q 2 )V s (k, q)O s (q, k') + (k' 2 - q 2 )O s (k, q)V s (q, k')} 
Jo 



(14) 



and decoupling is not manifest. Further, it is not clear if A provides a measure of decoupling 
in the case of a general operator. So we turn to visualizations of the operator matrix 
elements for guidance. As noted earlier, in practice we do not solve the flow equation for 
operators, but apply the unitary transformation Eq. Q at s by first solving the initial and 
final Hamiltonians for the eigenvectors and constructing the matrix 

Usfakj) = X)<W«( S )> <^(°)l fc ;> > ( 15 ) 

a 

where {k{\ is the discrete momentum mesh. 
B. Momentum Distribution 



We begin our examination of operator evolution with perhaps the simplest nontrivial 

i i 



example: the momentum occupation operator a\a in the center-of-mass frame (in this frame 
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FIG. 2. The momentum distribution in the deuteron as given by the expectation value of the bare 
number operator a q a q for the (a) Argonne v\% (AV18) pi] and (b) N 3 LO (500 MeV) [E] potentials, 
both evolved and unevolved. 



- for the A=2 particle space — q is the magnitude of both the relative momentum and 
the single-particle momentum). By varying the momentum q we can gain insight into how 
the SRG evolution behaves for initial operators dominated by either high or low momenta. 
In Fig. [2] the plot of the momentum distribution is reproduced from the expectation value 
(0d| a q a q \ipd) m the deuteron for the two initial potentials of Fig. |TJ The solid line is the 
result when the unevolved wave function (i.e., the wave function derived from the unevolved 
potential) is used with the unevolved a q a q operator. This sets the baseline for evaluating 
the effects of the SRG for each potential. When one uses the evolved wave function with the 
evolved operator, U s a q a q Ul, the lines are indistinguishable. The dashed and dot-dashed lines 
are calculated using the unevolved operator with the evolved wave function at A = 2.0 fm" 1 
and A = 1.5 fm -1 , respectively. These curves quantify the effect of not consistently evolving 
operators and also give the momentum dependence of the evolved wave functions. As can 
be seen, the high momentum components of the wave functions are significantly suppressed, 
as is consistent with the decoupling seen in the potential. 

An analogous visual representation to that used for the potential allows us to analyze the 
RG flow and properties of these objects. Figures [3] and [4] illustrate the flow of the momentum 
occupation operator consistent with the renormalization of the N 3 LO 500 MeV potential. 
We present the SRG results using only one potential, but qualitatively similar results will be 
obtained using any nuclear potential of interest. Each of the figures shows three sequences: 
the initial operator matrix elements (k\ a q a \k') evolved to four different A = s" 1 / 4 values 
and the integrand of (^d(s)| (a^ q a ) s {^(s)) with first linear and then logarithmic scales. 
The operators shown correspond to those used to calculate the momentum distribution at 
q = 0.34 fm -1 and q = 3.02 fm -1 (marked by the dotted lines in Fig. |2^b)). It is apparent 
that the unevolved operator is simply a delta function in momentum space (some minor 
evolution can already be seen at A = 6fm _1 because the scale must be magnified to view 
the evolution at lower values of A). 
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FIG. 3. (Color online) (a) SRG evolution of the operator (k\ a\a q \k') for q = 0.34 fm -1 in the 3 Si 
partial wave from A = 6fm _1 to A = 1.5 fm -1 , with the N 3 LO (500 MeV) [17] initial potential, (b) 
Integrand of (ipd(s)\ {a\a q ) s \4>d(s)) with linear (top) and logarithmic magnitude (bottom) scales. 



Consider the operator (top row) sequences first. For both q values, the evolution begins 
along the momentum axes around the delta function where we see strength developing 
that was not present in the original operator. This behavior can be understood from the 
momentum basis SRG [Eq. (14)] and the features of the corresponding potential evolution. 
Holding k fixed, one can see that only the second term in the integral initially picks up 
strength along the axis that passes through the delta function, because the operator is zero 
everywhere else, and vice versa holding k' fixed. For the operator at high q we see that 
it develops more and more strength at low momentum as it evolves. The need for this 
additional strength is particularly evident because of the decoupling via evolution of the 
potential and the consequent suppression of high-momentum components in the deuteron 
wave function, as seen in Fig.[2|b). As a result, the operator must pick up additional strength 
for the expectation values to remain unchanged. This strength can appear in two ways: (i) 
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FIG. 4. (Color online) Same as Fig. K3lbut for q = 3.02 fm 1 . 



It can come in at low momentum, as we see here; or, (ii) the operator can gain strength 
only at high momenta, where the operator would have to become pathologically large. If 
the second case were to occur, practical calculations with the SRG in a reduced basis would 
not be possible. This is found empirically to not be the case, and we can be confident such 
pathologies will not occur based on more general arguments discussed in what follows. 

The operator at low q, however, picks up some strength at larger values of momentum 
than present in the initial operator. This is also needed to compensate for the suppression of 
low momentum dependence in high-energy eigenstates to maintain their expectation values. 
One should note that the operator display scale used here can be a bit deceptive, in that it 
has been amplified to make the qualitative features of the evolution more apparent. Most 
of the evolution does, in fact, remain at low momentum for deuteron expectation values. 

We show the occupation operator as an integrand given by (-0d| a^ q a q \ipct) in Figs. p[b) and 
W[h). The expectation value filters the general operator by weighting its matrix elements 
with the deuteron wave function. Now we can see a clean RG flow in the strength for 
both operators. The integrand of the operator at q — 3.02 fm -1 begins as a sharp spike, 
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q [fm _1 ] 

FIG. 5. Decoupling in operator matrix elements is tested by calculating the momentum distribution 
in the deuteron after evolving the AV18 potential to several different A and then truncating the 
Hamiltonian and evolved occupation operators (i.e., set them to zero above A = 2.5 fm -1 ). 

corresponding to the original operator, but then flows out along the momentum axes to lower 
momentum. By the time the integrand reaches lower values of A in the evolution nearly all 
of the strength in the expectation value is in the low-momentum region. The original spike 
disappears as the wave function dependence at high momentum falls off. 

As for the operator at q = 0.34 fm _1 , the strength does begin to flow out to some extent 
but remains almost entirely in the low-momentum region. Once again, the display scale 
has overemphasized the extent of the evolution in the values of the integrand. The spike 
that remains at A = 1.5 fm -1 actually contains ~ 96% of the full expectation value. Owing 
to the possibility of misinterpreting these plots on a linear display scale, we also include 
the same plots with logarithmic display scale. These pictures show only the magnitude of 
the integrands, but display nearly the full range of their values. Now it is conclusive that 
the strength of the high momentum operator flows to low momentum, and the strength of 
the low-momentum operator remains at low momentum for a low-energy state. We see this 
pattern repeat in the calculation of other operators in the next section. 

Despite the apparent changes in the integrands as they evolve, it is important to note 
that the sum of all the points (the expectation value) remains unchanged because of the 
unitarity of the SRG transformation. The momentum distribution calculations shown in 
Fig. [2] were performed in the full momentum space of the original potential. Decoupling of 
the potential allows us to truncate the model space, thereby making numerical simulations 
more feasible, while at the same time allowing us to calculate the correct binding energies. If 
the calculation of other expectation values must be performed in the full model space, then 
the benefits of the SRG would be lost. However, the redistribution of strength implies that 
we have a form of decoupling for the operator. A critical test is to verify that decoupling is 
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q [frrf 1 ] 

FIG. 6. Decoupling in operator matrix elements is tested by calculating the momentum distribution 
in the deuteron after evolving the N 3 LO potential to several different A and then truncating the 
Hamiltonian and evolved occupation operators (i.e., set them to zero above A = 2.0 fm -1 ). 



maintained in the calculation of operator expectation values. 

This check is shown in Figs. [5] and [6] for the momentum distributions of AV18 and N 3 LO 
potentials. To perform the check, we evolved the original Hamiltonians and operators in 
the full momentum space to various A then truncated the model space to A. The deuteron 
wave function derived in this truncated space is used to calculate the expectation value of 
the evolved number operators to produce the momentum distributions shown here. The 
figures show that when the SRG evolution A > A, the curves deviate significantly from 
that produced in the full space (because the wave function is distorted). However, once the 
operators are evolved to A below the truncation at A, the expectation values are reproduced 
for all values of momenta, even in the region outside of the new model space. Thus decoupling 
is successful and A provides a rough guide as to where this decoupling occurs. We see that 
this is also the case for other operator matrix elements of interest, as well as understand 
further how this comes about, in what follows. 



C. General Analysis 

The plots of the deuteron integrands show that no pathologies appear at high momen- 
tum in the evolved operators and verify that decoupling can be successful when calculating 
expectation values of this operator. The plots even indicate where the model space can be 
truncated — this is simply where the integrand strength becomes negligible in the loga- 
rithmic plot. Here we develop a more general understanding of operator evolution to build 
confidence that pathologies will not occur in other operators. 
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Consider the representation of a generic operator in terms of the energy eigenstates, 

d s = J^o ij \Us))(^)\ , (16) 

where 

O tJ = (Us)\O s \^(s)) . (17) 

It is important to remember that these matrix elements are invariant under SRG trans- 
formations, so Oij does not depend on s. Thus, the momentum space behavior of evolved 
operators is given, in turn, by the momentum space behavior of the evolved eigenstates - 
specifically the sum of their outer products weighted by O^, 

O s (k,k') = J^Oij (k\^{s)) (^{s)\k') . (18) 

The behavior of these eigenstates is well under control. As we have seen from the momen- 
tum distribution of the deuteron, low-energy bound-state wave functions are suppressed at 
high momentum. The rest of the (positive) eigenstates are effectively smeared out delta 
functions (normalized in a finite basis), which because of decoupling in the Hamiltonian 
become increasingly narrow peaks with the evolution in s. So, the evolution will not become 
pathological unless the unevolved operator is already pathological (i.e., only if some are 
unnaturally large). 

If we now consider this operator in the deuteron eigenstate \ipd(s)) = \ipi(s)), we find that 
only the On matrix element of the operator is projected out and the momentum dependence 
is given by the outer product of the deuteron wave function. Specifically, the only nonzero 
part of the operator will be formally given by 

O s (k,k') ^ 11 (k\i; 1 (s)) (Ms)W) • (19) 

However, if we would like to reconsider the issue of decoupling and a finite model space trun- 
cation, we find that we are restricted by the potential breaking of eigenstate orthogonality 
in the truncated space; that is, the extent to which 

- [ A k 2 dk (^(s)\k) (k^jis)) ^ , for^j. (20) 

^ JO 

The logarithmic integrand plots show us that this problem is negligible if one of the states 
is the deuteron. Furthermore, SRG-driven decoupling of states well separated in energy 
will make them increasingly orthogonal in the truncated space. Thus, the momentum-space 
evolution of operators in a specific basis can be brought under control; we will consider the 
case of induced many-body components of the operators in a sequel to this article [10J. 

III. OTHER OPERATORS 

A. Long-distance Operators: (r 2 ), {Qd)i y ) 

We begin our presentation of additional operators with the evolution of operators for three 
paradigmatic expectation values in the deuteron: the rms radius, the quadrupole moment, 
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and r 1 . These operators all act on relatively long distance scales. At leading order they 
are naturally defined in coordinate space so that deuteron expectation values can be written 
as [IE] 

.1/2 

drr 2 (u{r) 2 + w(r) 2 ) 



(r d ) 

(Qd) = 



1 

" 2 
1 

20 



drr 2 w(r) (^\/8u(r) — w(r)j 



and 



dr 



w(r) 2 ) 



(21) 
(22) 

(23) 



where u and w are the 3 Si and 3 Di deuteron radial wave functions. However, we will continue 
to analyze the operators in momentum space. To avoid numerical instabilities associated 
with putting the derivatives in the momentum-space expressions on a mesh (see Ref. [18]). 
we extract the coordinate-space operators from Eqs. (21)-(23) and transform to the partial 



wave momentum basis. For example, from Eq. (21) we can see that the r operator is given 



by the diagonal matrix (discretized in coordinate space) 

(r| r 2 \r') = r 2 5{r — r') 
in the 3 Si and 3 Di channels. The only transformations needed are given by 

r 2 

and 



<r|fc; 3 Si) = J -rk 2 j (kr) 



(rlA^Dx) 



r k 2 jz(kr) 



(24) 



(25) 



(26) 



where the jj's are spherical Bessel functions (additional partial waves are needed for states 
other than the deuteron, of course). The transformations are represented asnxm matrices 
(where n and m are the sizes of the coordinate and momentum space meshes, respectively) 
and applied to both sides of the n x n coordinate-space operator matrix to produce an 
m x m matrix in momentum space for each partial wave. Then, to evolve the operator in 
momentum space, we apply the unitary transformation O s = U s O Uj. 
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FIG. 7. (Color online) Operator evolution of (k\ r 2 \k') in the 3 Si partial wave from A = 6fm 1 to 
A = 1.5 far 1 using the N 3 LO (500 MeV) potential. 
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FIG. 8. (Color online) Integrand given by {ipd\ r<1 iV'd) an d evolved from A = 6fm _1 to A = 1.5 fm -1 
using the N3LO 500 MeV potential with a linear color scale (top) and a logarithmic scale of the 
magnitude (bottom). Notice the difference in momentum scales. 



A 


6.0 fm -1 


3.0 fm -1 


2.0 fm -1 


1.5 fm -1 


Qd 


0.275 fm 2 


0.274 fm 2 


0.269 fm 2 


0.260 fm 2 



TABLE I. Expectation value of the quadrupole moment given by the unevolved operator with the 
evolved deuteron wave function. 



The SRG evolution sequence for r 2 shown in Fig. [7] is again a picture of just the operator 
in the momentum basis. Note the restricted momentum scale of the plot, as well as the 
magnitude of the operator display scale. As a long distance operator, the strength is highly 
concentrated at low momentum. It is evident that very little renormalization occurs at the 
lowest values of momentum, which is consistent with the findings for the number operator. 
Looking at the deuteron integrand of the r 2 operator in Fig. [8j the linear display scale 
plot shows the same behavior. The log-scale shows the entirety of the contributions to the 
integral. At low momentum, the strength is orders of magnitude greater than elsewhere. 
While the pattern of the contribution changes very little with A in the linearly scaled plot, the 
log-scale plot clearly shows the characteristic exponential suppression with decreasing A due 
to the decoupling of the potential and consequent reduction of high-momentum components 
in the deuteron wave function. 

The quadrupole moment operator is also a long-distance operator and we find its evo- 
lution shares most of the same characteristics found for the evolution of the r 2 operator 
(note that it also picks up strength in the 3 Si- 3 Di channel). The deuteron integrand is 
shown in Fig. [9j There is little change in the actual evolution of the operator, while high- 
momentum contributions become exponentially suppressed. The lack of evolution can be 
quantified by calculating the expectation value of the unevolved operator with the evolved 
wave function. Results for various A are given in Table [TJ Note that we are using the basic, 
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FIG. 9. (Color online) Integrand of the quadrupole moment expectation value in the deuteron and 
evolved from A = 6fm _1 to A = 1.5 fm -1 using the N3LO 500 MeV potential with a linear color 
scale (top) and a logarithmic scale of the magnitude (bottom). Notice the difference in momentum 
scales. 



one-body quadrupole moment operator without two-body or other higher order renormal- 
ized corrections. The "true" value for this potential and operator is Qd ~ 0.275 fm 2 while 
the experimental value of this quantity is Q exp t. ~ 0.285 fm 2 . Thus the induced two-body 
contribution is the same order as omitted two-body contributions to the initial operator. As 
such, the SRG has not led to any changes larger than one would expect from including the 
fully renormalized operator from the EFT. 

Finally, in Fig. 10 we show the evolution in the deuteron of the r~ l operator, which has 
larger contributions at short range than the previous examples. Consequently, this operator 
is more spread out in momentum space. Yet, we see the same general behavior with respect 
to renormalization and the suppression at large momenta without, as usual, any changes in 
the expectation value. 



B. Deuteron form factors: Gq, Gq, and Gm- 



We now turn to the SRG evolution of electromagnetic operators that determine the 
deuteron charge, quadrupole, and magnetic form factors (i.e., Gq, Gq, and Gm respec- 
tively) [T9H23] . We restrict our discussion to deuteron expectation values that have been 
derived consistently with chiral EFT at leading order in coordinate space |23j. These are 
given by 

G C (Q 2 ) = G ( £\Q 2 ) [dr[u 2 (r)+w 2 (r)]j (\q\r/2) , (27) 



G Q (Q' 2 ) = G ( £\Q 



6^2 

'IF 



dr 



u(r)w(r) + w 2 (r 



h(\q\r/2) 



(28) 



15 




0.5 

I 

* 1 



k' (fm" 1 ) k' (fm 1 ) k' (fm 1 ) k' (fm" 1 ) 

0.5 1 0.5 1 0.5 1 0.5 1 



V 


V 


v 




X =6.0 fm" 1 


X =3.0 fm" 1 


X =2.0 fm" 1 


A=1.5fm" 1 



u 





1 

il 



k' (fm" 1 ) k' (fm 1 ) k' (fm 1 ) k' (fm" 1 ) 

01234012340123401234 



FIG. 10. (Color online) Integrand given by (ipd\ r ~ 1 IV'rf) an d evolved from A = Qfm" 1 to A = 
1.5 fm -1 using the N3LO 500 MeV potential with a linear color scale (top) and a logarithmic scale 
of the magnitude (bottom). Notice the difference in momentum scales. 




G ( v\Q 2 )l f drw 2 (r) [j (\q\r/2)+j 2 (\q\r/2)] 
+ G$(Q 2 )2 Jdru 2 (r) Jo (\q\r/2) 
+ G$(Q 2 ){V2 fdru(r)w(r)j 2 (\q\r/2) 

drw 2 (r)[ j0 (\q\r/2) - j2 (\q\r/2)]} 



(29) 



where Q 2 = \q\ 2 , and G^>(Q 2 ) and G\°>(Q 2 ) are the sinele-nucleon isoscalar electric and 
magnetic form factors obtained from the parametrization given in Ref. [21]. From these 
coordinate-space expressions, we can apply the same procedure used earlier to extract the 



operators, then use Eqs. (25) and (26) to convert to momentum space and transform via 
the SRG unitary transformation. One should note that starting from a coordinate-space 
operator is by no means essential; it is simply a numerical convenience in this case. 

The expectation values as functions of q = x/Gp for each of these operators is presented in 



Fig. 11 The solid line has been calculated using the unevolved potential with the unevolved 
operator, which again serves as our reference value. The starred points are calculated using 
the evolved wave function with the evolved operator (both at A = 1.5 fm -1 ). As advertised, 
they lie precisely on top of the solid line for all values of q, up to small numerical errors. The 
dot-dashed line is calculated using the unevolved operator with the evolved wave function 
as an indication of the effect of renormalization on the expectation value. We see noticable 
deviation above q ~ A; however, from the magnitude of the suppression seen in the wave 
function at high momentum one might have expected the curve to drop much faster with 
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FIG. 11. Deuteron form factors Gq, Gq, and Gm using the isoscalar electric form factor 
parametrization from Ref. [24J. M is the nucleon mass and is the mass of the deuteron. 
The wave function is derived from the NNLO 550/600 MeV potential and the evolution is run to 
A = 1.5 fm -1 ESI- 



respect to q. However, the form factor operators probe momenta in the deuteron center-of- 
mass frame whereas q is specified in the laboratory frame. Thus, the operators are probing 
the wave function largely at ^q, which again brings the calculations in line with our SRG 
expectations. 

The basic features of the SRG evolution of all three operators are qualitatively very 
similar, so we present the visual matrix representation of the magnetic form factor at high 
and low q as a representative example in Figs. 12 and 13 This form factor picks up strength 
in all deuteron channels. The high-momentum form factor has a much greater diffusion of 
strength at high momentum than seen in any of the static properties explored earlier. Yet 
it is apparent that the strength in the operator flows to low momentum in this case also. 
In contrast, the low-momentum operator exhibits very little renormalization, as we have 
come to expect. For both cases, we can see from the logarithmic display scale plots that the 
momentum dependence of the form factors is virtually eliminated at large momenta without 



17 



affecting the outcome of the computation. 
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FIG. 12. (Color online) Integrand of Gm at q = 6.90 fm -1 evolved from A = 6fm _1 to A = 1.5 fm -1 
using the NNLO 550/600 MeV potential with a linear scale (top) and a logarithmic scale of the 
magnitude (bottom). 



This lack of dependence on high momenta in evaluating expectation values is particularly 
significant for the practical application of the SRG in calculations of low-energy few- or 
many-body systems. Had this renormalization led to singular behavior in the operators at 
high momentum, the effects due to the suppression of the wave function would have been 
negated and led to wildly erroneous results in the evaluation of observables in a reduced 



model space. The arguments in Sec. II C explain why this will generally be the case. 



IV. VARIATIONAL CALCULATIONS 

We have argued previously that SRG-evolved interactions and the resulting wave func- 
tions become "simpler." Variational calculations of the ground-state wave function can pro- 
vide a test (and additional meaning) of the extent to which this is true in practical ap- 
plications. The decoupling of high- and low-momentum states caused by SRG evolution 
means that the resulting wave functions are much less correlated than the original wave 
functions. Consequently, one expects that variational calculations of the evolved wave func- 
tions should be effective with a much simpler ansatz than one would normally require for 
the corresponding unevolved interaction. In conjunction with this, one might be concerned 
that a delicate interplay of the evolved operators and wave functions would be necessary 
to preserve matrix elements. However, it turns out that the evolved operators are not only 
equally good, but actually superior to the original operators in this respect. Rather than 
high-momentum operators picking up small pieces of the wave function (which could never 
be reproduced by a simple variational calculation), we get a smooth sum over where the 
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FIG. 13. (Color online) Integrand of Gm at q = 0.34 fm -1 evolved from A = 6fm _1 to A = 1.5 fm -1 
using the NNLO 550/600 MeV potential with a linear scale (top) and a logarithmic scale of the 
magnitude (bottom). 
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wave function is large and easily approximated. We illustrate these points here by choos- 
ing a simple variational ansatz and looking at how convergence improves for SRG-evolved 
operator expectation values. 

Variational calculations have been performed on low-momentum potentials (specifically 
"Mow*:," derived using an alternative RG formulation) in the past to demonstrate a sig- 
nificantly improved convergence of the binding energies [26j |27]. We choose to adapt a 
simple ansatz for the deuteron used in those calculations to make our point for the SRG. In 
particular, we take the (unnormalized) 3 Si and 3 -Di partial waves to be 

u ( k ) = 777 77777 97 e ^ > ( 30 ) 

v ' (k 2 + 7 2 )(/c 2 + /i 2 ) v ' 

ak 2 f k2 \ 2 

= jjF^FW^f 6 ' (31) 

where 7, fi, u, and a are variational parameters. The exponential factors are chosen to 
match the asymptotic suppression of the wave function resulting from the decoupling of the 



interaction according to Eq. (12). The energy is minimized with respect to the variational 
parameters at various A for the three different potentials used in this article. The binding- 
energy results are shown in Fig. 14 Without evolution, the AV18 and N 3 LO trial wave 



functions are not even bound, and the NNLO wave function accounts for less than half 
the binding energy. With evolution, the AV18 and N3LO wave functions begin to bind at 
A ~ 7fm _1 and A ~ 4.5 fm -1 respectively and when the evolution is taken further, the trial 
wave function is able to reproduce the exact binding energies to within 1 keV. 

Examining the matrix elements of operators which initially have strength concentrated 
over a range of different momenta — such as the occupation operators with respect to 
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FIG. 14. Deviation from of the best variational energy as a function of SRG decoupling 



parameter for the wave function ansatz of Eq. (31). Ed is the deuteron binding energy for each 



interaction derived via a full eigenvalue calculation of the Hamiltonian. 
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FIG. 15. The momentum distribution in the deuteron as given by the expectation value of 
the evolved occupation operator Ua\a q U^ using the variational wave functions derived from the 
Salpeter ansatz and the AV18 potential evolved to A = 6.0 fm -1 , 4.0 fm -1 and 1.5 fm~ . The direct 
calculation is from a full eigenvector solution of the Hamiltonian. 



momenta q — provides a stricter test of the variational solution to the evolved wave functions 
and the sensitivity of evolved operators to them. The initial AV18 potential has particularly 
strong correlations at high momenta. If we look at the evolved occupation operators in 



Fig. 15, we see three curves: one with a wave function near the binding threshold, one at 



about half the binding energy, and one that is well converged with respect to the binding 
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FIG. 16. Deuteron form factor Gm using the isoscalar electric form factor parametrization from 
Ref. [23]. M is the nucleon mass and is the mass of the deuteron. The variational wave 
functions are derived from the Salpeter ansatz and the NNLO 550/600 MeV potential evolved to 
A = 5.0 fm -1 and A = 1.2 fm . The operators are evolved consistently. The direct calculation is 
from a full eigenvector solution of the Hamiltonian. 



energy (evolved to A = 6.0 fm -1 , 4.0 fm -1 and 1.5 fm -1 respectively). Near the binding 
threshold, the momentum distribution is reproduced rather poorly, but at smaller A the 
curve improves, and once the binding energies are converged the operator expectation values 
are also converged to approximately 1% or better. 

The same pattern holds for other operators and interactions; that is, the operator matrix 
elements are not sensitive to the fine details of the evolved wave function. The magnetic 



form factor of the deuteron using the NNLO potential, for example, is shown in Fig. 16 Not 



only does this operator pick up strength in both partial wave states of the deuteron, but also 
their coupling. Again, for a variational wave function at half binding energy (A = 5.0 fm" 1 ) 
the matrix elements deviate significantly from the direct, no n- variational calculation, but 
when the binding energy is converged (at A = 1.2 fm^ 1 ), the form factor expectation values 
are reproduced to better than 1%. 



V. OPERATOR FACTORIZATION 



In this section, we consider in more detail the expectation value in a low-energy bound 
state of operators that initially have strength only at high momentum. The momentum 
distribution of the deuteron at large q is our prototype. The momentum distribution for the 
initial potentials in Fig. [2] show structure at high momenta because of the short-range repul- 
sion in the potential (particularly for AV18). When evolved to low momentum, this structure 
disappears and the deuteron wave function will select the low-momentum part of the evolved 
operator. But this evolved operator must still reflect the external high-momentum scale. We 
can anticipate simplifications by exploiting this separation of scales provided by the SRG; 
in particular, we expect a factorization of the evolved operator based on operator-product- 
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expansion (OPE) arguments applied to nonrelativistic effective theories [28j 129 



A. Numerical Verification of Factorization 
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FIG. 17. Numerical tests of factorization of the unitary transformation U x (k,q) by plotting the 



ratio in Eqs. (32) and (33) as a function of q for fixed and several values of ki- Plateaus in q 
indicate factorization. The unitary transformations are generated from the Argonne vig (AV18) [16] 
potential evolved to A = 2 fm -1 in the (a) 1 So and (b) *Pi partial waves. The shaded region marks 
q < A. 



Previous calculations of the deuteron momentum distribution suggested that the unitary 
evolution operator, U x (k, q)^factorizes into a function of k times a function of q, U\(k, q) — > 
K\(k)Q\(q), for k < A and q ^> A [13J. To numerically test for factorization in the unitary 



transformation, we use transformations generated via Eq. (15) and the evolution of NN 
potentials, and consider the ratio 



Ux{k u q) , K x (k t )Q x (q) 
U x (k ,q) ~^ K x (k )Q x (q) 



(32) 



holding hi and k constant with k A. If there is factorization, the q dependence should 
cancel; that is, for k < A and g > A we should find 



U\{ki, q) K x (h 



U x (k ,q) K x (k ) 



(33) 



In Figs. 17 and 18 we plot the ratio in Eq. (32) versus q for representative cases. The 



signature of factorization is a plateau in q. The shaded regions are where q < A. In all 



2 Because A is an important momentum scale in our factorization discussion, we use the notation U x with 
A = s^ 1 / 4 rather than U s in this section. 
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FIG. 18. Same as Fig. 17 but for the 3 Si partial wave and A of (a) 1.5 fm 1 and (b) 3fm 1 . 



cases, there is no factorization in this region, consistent with the requirement that q ^> 
A. In the unshaded region we see definite plateaus for q > A as long as fcj < A, with 
diminishing prominence as ki increases (they disappear for hi > A). Thus we have at least 
a qualitative verification of factorization. Note that Fig. [l8](b) shows that for larger A the 
clean factorization breaks down (as well as restricting the applicable domain). 

The singular value decomposition (SVD) can be used as a tool to quantitatively analyze 
the extent to which U\ factorizes. The SVD of a matrix M can be expressed in general as 
an outer product expansion 



M 



(34) 



i=l 



where r is the rank of the matrix and the di are the singular values (in order of decreasing 
value). The idea is that if the first singular value, d\, is sufficiently large compared to the 
others, the first term dominates and we have a factorized approximation. We can apply 
this to U\ in the region where high and low momentum couple. Thus, the vector ui would 



correspond to the low momentum function K\{k) from Eq. (32) and Vj to Q\(q). If valid, 
one can calculate the factorized operator using the unitary transformation obtained directly 
from the SVD. Moreover, the expansion provided by the SVD allows us to make systematic 
corrections to the factorized unitary transformation and the operators evolved with it. 

To test if such an expansion can be used, the first few singular values have been calculated 
in Table [TT] for the q > A and k < A region of the SRG unitary transformations for several 
different potentials, each evolved to A = 2fm~ 1 . That is, the SVD is applied to the matrix 
obtained when elements of U\(k,q) with k > A and q < A are set to zero; in practice a 
cutoff A « 2.5 fm -1 is used to ensure that we are in the region where off-diagonal coupling 
is strongly suppressed everywhere in the Hamiltonian. The dominance of d\ in each case is 
promising. 

To test if a truncated outer product sum is a good approximation to the contribution from 
k < A, q > A, we consider the errors in some representative expectation values in Table [TTT 
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3 Si- 3 D! 


Potential 


di 


d 2 


ds 


di 


d 2 


dz 


AV18 


0.763 


0.033 


0.007 


0.671 


0.015 


0.008 


N3L0 500 MeV 


1.423 


0.221 


0.015 


1.873 


0.225 


0.044 


N3LO 550/600 MeV 


3.074 


0.380 


0.061 


4.195 


0.587 


0.089 



TABLE II. Singular values of the unitary transformation U(k,q) for q > A and k < A (see 
discussion in text; units in fm -1 ) corresponding to the given potentials at A = 2 fin -1 in the 1 So 
partial wave and 3 Si- 3 Di coupled channel. 

for several levels of truncation. The "zeroth-order" contribution is from the matrix where 
k > A and q > A is set to zero (this is denoted by SVD=0 in the table). The first-order 
(SVD=1 in the table) contribution uses the same matrix plus the approximation of U\(k, q) 
for k < A and q > A given by the first outer product in the SVD expansion. The second 
order approximation uses two outer products, etc. The occupation and charge form factor 
operators shown here are chosen to illustrate the effects of the factorized approximation at 
various momenta. Additionally, the initial occupation operator has no off-diagonal strength, 
whereas the initial charge form-factor operators have relatively substantial off-diagonal con- 
tributions at large values of q; this is significant for the applicability of factorization to an 
operator, as we see below. So, what we find at low momentum for the occupation operator is 
that it is essentially exact, up to errors resulting from decoupling and truncation of the wave 
function, and it is the same with or without the SVD approximation. This error increases, 
as expected, for larger A. Because Gq is more diffuse initially, we see a small improvement 
even at small momenta when using the SVD approximations. 

For the occupation operator at high momenta (well above the cutoff), the error is 100% 
without an approximation to U\(k, q) because there is a hard cutoff and the initial operator 
is localized in the upper region of momentum space. However, with just one term in the 
SVD expansion we recover that expectation value for A = 1.5 fm -1 to better than 1%, and 
the situation improves with additional terms in the expansion. At A = 3.0 fm -1 decoupling is 
evidently not sufficient for this approximation to work well. At very large values of momenta 
(e.g., q ~ 6.9 fm^ 1 ) the charge form factor shows improvement with the SVD approximation, 
but because this operator has significant off-diagonal strength, the improvement is not as 
pronounced. At a value of q ~ 3.0 fm _1 it is evident that the SVD still improves the relative 
error. However, recall that the strength in the form factors is larger around ~ \q. 

B. Connection to the Operator Product Expansion 

The OPE was developed for the evaluation of singular products of local field operators at 
small separation. In our case, where such operators are treated as matrices and we typically 
work in momentum representation, the focus becomes low-momentum matrix elements of 
a product in which high-momentum states dominate the intermediate sum. This leads us 
directly to consider low- to high-momentum matrix elements of SRG-evolved operators, and 
a generic analysis is then based on the study of U\(k, q) for k < A and q ^> A. 

The utility of the OPE rests on factorization; short-distance details decouple from long- 
distance dynamics. Factorization enables one, for example, to separate the momentum and 
distance scales in hard-scattering processes in terms of perturbative QCD and parton dis- 
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SRG A = 1.5 fin -1 


SRG A = 3.0 fnr 1 


Operator 


SVD 


g = 0.34 


q = 3.01 


q = 6.90 


SVD 


q = 0.34 


q = 3.01 


q = 6.90 


(4 °>q) 





7.61xl0~ 7 


1.00 







1.06xl0~ 3 


1.00 




with N 3 LO 


1 


7.61 xlO" 7 


4.28xl0" 3 




1 


1.06xl0" 3 


6.36X10" 1 






2 


7.61 xlO" 7 


4.79 xlO" 4 




2 


1.06xl0" 3 


6.35X10" 1 




G C (q) 





6.90xl0" 4 


5.01xl0~ 3 


8.93X10" 1 





4.10xl0" 4 


3.36xl0" 3 


8.92xl0 _1 


with NNLO 


1 


1.28xl0" 7 


8.90xl0" 5 


4.06xl0" 2 


1 


1.63xl0" 4 


2.66xl0" 4 


4.00X10" 1 




2 


1.04xl0" 6 


2.10xl0~ 5 


4.18xl0~ 2 


2 


1.63xl0" 4 


3.04xl0" 4 


4.09X10" 1 



TABLE III. Relative error of evolved operator matrix elements calculated using the SVD to 
factorize U\(k, q) in the region where k < A and q > A (see discussion in text; units in fm _1 ). 



tribution functions. In our case, factorization is the direct result of decoupling. It provides 
tools that let us parametrize the high-momentum components of operators that would nor- 
mally require degrees of freedom we do not retain. We can, for example, build effective 
few-body operators containing state-independent functions of high momenta that can be 
measured directly in few-body experiments. These operators can then be employed to make 
predictions for A-body systems. 

Consider a generic operator, 0\ = U\OU\, and employ the spectral representation for 

U x : 

U x (k,q) = £<A#a(A))<Va(oo)|g). (35) 

a 

The OPE deals with cases in which the unevolved operator is dominated by high momenta 
(e.g., a) q a q with large q is the simplest paradigm) and we focus on k < A and q ^> A. For 
k < A we exploit the fact that low-momentum components of high-energy eigenstates of H\ 
are exponentially suppressed because of decoupling. As a result the sum is dominated by 
low-energy states, 

U x (k,q)K (k\M*))(M°°)\Q)- ( 36 ) 

E a -€.\ 2 

Once the sum is restricted we can turn our focus to approximating the high-momentum 
components of the unevolved low-energy states. This analysis is closely related to Lepage's 
discussion of the OPE analysis of wave functions, which leads him to write for S-waves in 
position representation [28J: 

^truc(r) = 7(0 f d 3 r ^ cS 5 3 a (r) + r){r)a 2 J d 3 r ^ cS V 2 5 3 a (r) + 0{a A ) , (37) 

where the coefficient functions 7(r) and fj(r) are state-independent parametrizations of the 
short-distance physics, and a is approximately the distance of the ultraviolet cutoff. In 
this section, we outline how SRG factorization can be understood more generally (and 
analytically) in the context of the OPE for nonrelativistic Schrodinger problems by deriving 
an analogous equation in momentum space for the SRG-evolved wave function. 
To do so, we first define the projection operators 

Va= [ A dp\p)(p\ (38) 
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and 

Qa= I dq\q)(q\, (39) 



where A divides momentum space and dp = dp in the partial-wave momentum basis. 
This A is to be distinguished from A, which is the SRG evolution parameter and an approxi- 
mate measure of decoupling in the evolved potential. We use ip* to denote the eigenstates of 
the Hamiltonian ordered according to increasing energy E a and evolved to A via the SRG. 
H\ and V\ represent the corresponding SRG evolved Hamiltonian and potential. The initial, 
unevolved operators correspond to A = oo. 
From the unevolved Schrodinger equation 

ffoo |C> = E a , (40) 

we can write 

(VaH^Va VaH 00 Qa\ (VaK\ _ F (Vk^\ uu 

{QaH^Va QaH^Qa) \Qa^J a \QA^J ' 1 } 

and thus for the "Q" space we have 

Qa |C> = (E a - QaH^Qa^QaH^VaVa |C> 

= (E a - QaHooQaVQaVooVa , (42) 

where we have used ("Pa) 2 = Va, #oo = T + Voo, and QaTVa = 0. For low-energy states 
t/>£° such that \E a \ <C Min[|£ , Q_H-g|] (where Eqhq are the eigenvalues of QHQ), we can 
neglect the E a dependence. Also, assuming that the potential V^q' ,p) is slowly varying 
with respect to p compared to Vq°(p) m ^ ne region p < A and q' ^> A, we can use the 
expansion for S-waves 

A nA 

dpV oo (q / ,p)^(p)^V oo (q / ,p')\ pl=0 x / dp^{p) 



o 



d 2 



A 

a„2 



x / dpp'^{p) +■■■ (43) 

p'=0 



to first order, combined with the fact that the low-energy states will have momentum com- 
ponents peaked at small p, to write 

roc i* A -i 

<?l« « " / dq' / dp (q\ \q') VM, 0) ^(p) . (44) 

J A JO yA-noo^A 

Tests indicate that these assumptions are valid for realistic AW potentials. 

Further, we see empirically via Fig. [2] that Va |Vo°) ~ when A > A (this is 

consistent with our understanding that the SRG with G s = T re i renormalizes/suppresses 
only the short-distance components of the wave function for values of A considered here). 
Thus, setting A = A and defining 

poo -I 

7 A (g) ee - / dq' (q\ ^-—^ \q') V^q',0) , (45) 



we have 



^a(q)^i\q) I dpz(X)ij x a (p). (46) 
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So we see that the high-momentum components of low-energy eigenstates can be factorized 
into a state-independent function 7 A (g), which summarizes the short distance behavior of the 
wave function, and a coefficient (given by an integral over the renormalized wave function) 
that gives the contribution due to the long-distance structure of the state. Moreover, if we 
include higher-order corrections resulting from the expansion of jQdpV 00 (q',p)ilj'^ ) (p) about 
p = 0, we recover the analog to Lepage's OPE, Eq. (37), in momentum space for the short- 
distance structure of a wave function. It is given by 



<7>/i ^ l\q) l X dpZ(\)ip x a (p) + V \q) f X dpp 2 Z{\)4> x a {p) + -- 



where 7 A (g) is given previously and 



(47) 



(48) 



p=0 



Now, from the definition of the SRG unitary evolution operator in Eq. (15), in the region 
k < A and g> Awe can use the leading-order term of our OPE to write 



Ux(k,q) = J2( k \^) (CI?) 



\E a \<.\E QHQ \ 



Z / dpZ{X)^\p) 

= K\{k)Q\(q) , 



7 A (?) 



(49) 



where the sum is only over states in the "P" space thanks to decoupling. Thus, we can 
understand the factorization of U\ as a general consequence of our ability to factorize the 
high-momentum components of low-energy nuclear wave functions via an OPE plus decou- 
pling in the SRG-evolved Hamiltonian. Moreover, because i[>a(k) to good approximation has 



no support when k < X for a in the "Q" space, we can extend the a sum in Eq. (49) to the 
full space and apply closure to find 



Ux(k,q) 



Z(X) / dpJ2(k\i>a)(i>a\p) 



7 A (g) « Z(A) 7 A (g) • 



(50) 



Thus, to a first approximation, K\(k) is a constant factor. 

This approximate constancy implies that the ratios for the L = channels in Figs. |l7^a) 
and [l8^a) and (b) should tend to one in the factorization region, which is realized at the 
10-20% level for sufficiently low A. For L > 0, the generalization of Eq. (50) follows from 



modifying the Taylor series in Eq. (43) to account for Voo(q',p) oc p L for small p. Then 
Eqs. (45) and (46) are changed to 



i\q) 



1 d L 
QxHooQx dp L 



p=0 



and 



A 



«9)~7 A (9) / dpZ(X)p L ^ a (p) 



(51) 



(52) 
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With these changes, Eq. (50) for the factorization at leading approximation becomes 



U x (k,q) 



Z(X) / p L dpJ2(k\^)(^\p) 



7 A (g) « ^Z(A) 7 A (g) • 



(53) 



This approximation implies that the ratios in the factorization region should tend for L > 
to (ki/ko) L , which is seen at the same 10-20% level in Fig. 17 b). 

To gain insight into the implications of this factorization, we consider the expectation 
value of a q a q in a low-energy state, the deuteron. Because we know that strength in the 
evolved number operator expectation value decouples from high-momentum contributions 
in the deuteron, we can write 

<^ A | (a\a q ). |^ A ) = <^ A | U x (a\a q ) U\ |^ A > 



dk> 



dk' 



dq' dq" dk ^ (k')Ux(k' , q')8(q' - q)S(q" - q')U x (q", k)%j)^(k) 



(54) 



dk ^(k')Ux(k',q)Ux(q,k)^ x d (k). 



For a low-momentum operator, one with q < A, the expectation value thus depends only on 
the low-momentum details of the wave function (original and evolved). For q ^> A, however, 
we can make use of factorization and set U(k, q) — > Kx{k)Qx(q) to write 



dk' 



dk^ x d \k')Kx(k')[Qx(q)Qx(q)]Kxm X Ak) 



(55) 



from Eq. (54). Here we see that the expectation value of a high-momentum number operator 



is independent of the long-distance structure of the wave function. This is consistent with 
earlier calculations of the deuteron momentum distribution [13]. Again, as with decoupling 
in the potential, we appear to have a means by which long- and short- distance details can 
be separated for an operator evolved via the SRG. 

The generalization of this result is straightforward. Consider the expectation value of 
an arbitrary operator, 0(q',q), in a low-energy state, Vw- Because decoupling is valid for 
operator expectation values in a momentum basis (as we have seen via the expectation value 
integrand plots in Sees. [TT] and [TTT]) , we can write 



inw I u x ou x 



dk' dq' dq dk[^k')^Ux(k',q')0(q',q)Ux(q,k)^ ow (k) 



(56) 

We separate the integrals over the operator in the expectation value and apply factorization 
to set U(k, q) — > K\{k)Q\(q) in the region where k < A and q ^> A. If the unevolved operator 
has coupling between high and low momentum above the the factorization cut, then there is 
no great simplification. However, if the unevolved operator does not have coupling of high 
and low momentum above the factorization cut, factorization will allow us to separate out 
the high- and low-momentum structure of an operator into two contributions: 



dk' / dq 



dk lri ow (k')VUx(k', q')0(q', q)Ux(q, k)^ ow {k) 



(57) 
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and 



dk> dq> dq / dM^ w (^)] t ^A(A; / )[QA(g / )C > (? / ,g)gA(9)]^A(A;)V'i A owW- & 



This is analogous to what was found for the number operator. 

Thus, we see that the breakdown of contributions to the expectation value of a general 
operator is consistent with our interpretation of the SRG flow equations as a means by which 
one can achieve a separation of scales in the evaluation of nuclear few- and many-body prob- 
lems. We see explicitly here that the effects of a low-momentum probe of the ground state 
wave function depends (almost entirely) on the low-momentum details of the renormalized 
wave function. Likewise, the effect of a high-momentum probe is largely independent of the 
low-momentum structure. It is only for operators which probe the coupling of high and low 
momentum (long and short distance) details of the wave function that we must consider the 
full momentum space evolution of the operator. For the operators that have been considered 
in this paper, the latter is only true of the electromagnetic form factors at relatively high 
momenta (beyond the typical regime of interest for nuclear structure); for any operators 
which weakly couple high and low momentum, these terms can be neglected. 

To summarize, we can write the expectation value of an operator that has weak coupling 
between high and low momentum as 



<^ A ow | u x ou{ |^ A ow > « / "dk' tdk [^L(k')] j 

Jo Jo 



A p\ 

dq 1 I dq 



Uo Jo 

x U x (k',q')0(q',q)U x (q,k) + I QOQ K x (k')K x (k) 



< w (k) ,(59) 



where 



/*OG /*00 

I Q OQ= dq' dqQ x (q')0(q',q)Q x (q) . (60) 
J\ J\ 

By using factorization, we have seen that the expectation value breaks into a sum of two 
components: one which describes low-momentum structure, and a high-momentum compo- 
nent that factorizes into a piece depending on the low-momentum structure, and another 
piece that is a universal function of high momentum q. 



C. Interpreting high- momentum scaling behavior of the momentum distribution 

As discussed in Sec. [TJ results from large-momentum-transfer experiments such as (e, e'p) 
have been related [15J to how the tails of momentum distributions calculated for nuclei 
and nuclear matter using phenomenological potentials exhibit a scaling behavior at high 
momentum [30J. In particular, the momentum dependence of the distributions for nuclei 
ranging from the deuteron to oxygen, as well as nuclear matter, at high momentum is similar 
except for an overall nucleus-dependent scaling factor. One explanation for this is based on 
the dominance of two-body forces in the interaction and short-range correlations in the wave 
functions [T5] . How can we explain this feature in an SRG-evolved calculation, for which 
high-momentum components and short-range correlations are suppressed? 

Factorization provides a compelling alternative explanation. Because the a q a q operator 
for large q has no coupling to low momentum, the entire q dependence comes through the 



function Iqoq in Eqs. (59) and (60), which is independent of the low-momentum part. If 
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induced many-body contributions to the operator are relatively small and we neglect for the 
moment effects from embedding two-body operators into an A > 2 space, we conclude that 
for A > 2 the momentum distribution should be approximately the same for every A, with 
a scaling factor given by an A-dependent low-momentum integral over the low-energy wave 
functions. 




FIG. 19. The scaling of momentum distributions at high momenta in a ID model is tested by 
using the leading-order factorized approximation to the momentum occupation number operator 
to predict high-momentum scaling in A = 2,3,4 (symbols). The full momentum distributions for 
A = 2, 3, and 4 are shown with solid, dot-dashed, and dashed lines, respectively. 



We have tested this proposal in a ID model with an interaction that mimics features 
of the nuclear NN potential. The model and procedures used for the ID calculations are 
described in Ref. [12]. The full momentum distributions for two, three, and four particle 
systems in this model are shown with solid, dot-dashed, and dashed lines, respectively, in 



Fig. 19 The behavior at high momentum is analogous to the nuclear calculations [30J: 
The momentum dependence is similar for each system so that each curve differs only by a 
scaling factor. We then evolve the model interaction via the SRG to A = 2 and extract the 
unitary transformation. Only the operator from A = 2 is used; that is, any induced three- 
or four-body component is neglected. Details of the extraction and embedding of operators 
for A > 2, including boost corrections, are described in the sequel [10]. By using the first 
term in the SVD expansion to factorize U(k, q) in the region where k < A and q ^> A, we are 
able to reproduce to a large extent the momentum distributions at high momenta (shown 



with symbols in Fig. 19), and confirm our expectations regarding the scaling of the curves. 



These results are very promising and merit further investigation. 
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VI. SUMMARY 



In this article, we have examined the evolution of operators via the SRG with re- 
stricted application to the deuteron. We considered only the most commonly used generator 
G s = [T re i, H s ] with normal ordering in the vacuum. At this two-particle level it is easy to 
ensure that the transformations are unitary to high accuracy, so the invariance of matrix 
elements is assured. Thus our focus is instead on the nature of the evolved operators: Does a 
form of decoupling apply? Do operators become increasingly complicated as the wave func- 
tions become increasingly less correlated? How large are induced two-body contributions to 
various one-body operators? 

By considering the operator matrix elements in momentum representation both with 
and without deuteron wave functions included, we are able to follow the flow of strength. 
Because the transformations are unitary, the integrated value does not change with A, but 
the nature of the operator does. There is little evolution in long distance operators, whereas 
high-momentum operators must evolve significantly to compensate for suppression of high 
momenta in low-energy wave functions. In the end, one can see that the movement of the 
strength in the operator expectation values is given by the evolution of the eigenstates of 
the Hamiltonian itself. Moreover, we find that decoupling succeeds for operator expectation 
values in general, not just for the binding energies. 

The momentum distribution is particularly interesting because the evolution of high- 
momentum operators leads to their strength flowing completely to low momentum. Thus 
while the deuteron wave function has rapidly decreasing support at the high momentum, 
its matrix elements are preserved without pathologies in the transformed operators. In- 
deed, operator matrix elements are less sensitive to details, as evidenced by the improved 
effectiveness of estimates using variational wave functions. Decoupling for operators fol- 
lows as the contributions from higher energy/momentum basis states become unimportant, 
allowing truncation. This was explicitly illustrated in Figs. [5] and [6] The generality of 
these conclusions is evident by considering the eigenvector expansion of the SRG unitary 
transformations, which dictates the flow of strength. 

For low-momentum operators, which also includes the low-momentum part of one-body 
electromagnetic form factors, there is relatively little running and therefore only small in- 
duced two-body parts (which for A = 2 is simply the difference between the initial and the 
evolved result). In general, if the initial operator matrix elements pick up their strength 
predominantly at long distance, the operators will evolve only slightly until A is small. For 
electroweak operators, the real interest is in few- and many-body systems, where the simpler 
SRG-evolved wave functions are most advantageous. 

The practical pathway to few-body operators and applications to A > 2 is at present 
through the (Jacobi) harmonic-oscillator basis. The successful SRG evolution of three- 
body forces (and higher in model systems) in this basis is detailed in Refs. [El [12] • The 
corresponding challenge is to evolve operators for A = 2 and A = 3 and then embed them 
(including induced contributions) in higher- A spaces. Procedures for carrying this out, 
including the need to boost some operators, will be discussed in the sequel to this article 
[TO] . The calculation of transition matrix elements rather than ground-state expectation 
values requires separate consideration. 

Of particular interest will be the further study of SRG operator factorization, which 
occurs when there is a scale separation between the initial (unevolved) operator and the 
wave function momentum scale, which is limited by A. This factorization was shown to be a 
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natural consequence of applying the OPE to the unitary transformation. The extension to 
A > 2 was previewed in an application of factorization for momentum distributions of low- 
energy bound states, which provided an alternative interpretation to the commonly invoked 
role of short-range correlations. Work is in progress on the full realistic three-dimensional 
calculations. 

In closing, we reiterate that the favorable consequences of the SRG is a specific realization 
of the more general observation that the RG allows one to focus on the most relevant degrees 
of freedom in a physical problem [31 J. Thus, an evolution to low momentum for nuclear 
systems can be win-win not only for the Hamiltonian and wave functions but for operators 
as well. The SRG has some special advantages in practice because it uses operator flow 
equations that can be applied in any convenient basis with a variety of options for tailoring 
the flow. Furthermore, residual dependence on the flow parameter s or A becomes a powerful 
tool for assessing approximations. Different resolutions can lead to very different physical 
interpretations and intuition; the RG has the great advantage of being able to connect the 
different pictures. 
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